{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Back to the widget production problem\n",
    "Let's go back to the Opti 101/201 classic: Widget production and distribution. Use the code below as the base model for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*} \n",
    "{\\rm minimize} \\space &\\sum_{p,d}c_{p,d}x_{p,d}, \\forall p \\in P, d \\in D\\\\ \n",
    "\\sum_{p}&x_{p,d} \\ge n_d, \\forall d \\in D \\\\ \n",
    "\\sum_{d}&x_{p,d} \\le m_p, \\forall p \\in P \\\\ \n",
    "\\sum_{d}&x_{p,d} \\ge a*m_p, \\forall p \\in P\\\\ \n",
    "&x_{p,d} \\in \\{0\\} \\cup [30, u_p]\n",
    "\\end{align*}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    <strong>Note!</strong>\n",
    "    <p>Make sure to run the next two code cells to get a baseline for the original problem.</p>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in transportation cost data\n",
    "path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_1/'\n",
    "transp_cost = pd.read_csv(path + 'cost.csv')\n",
    "\n",
    "# get production and distribution locations from data frame\n",
    "production = list(transp_cost['production'].unique())\n",
    "distribution = list(transp_cost['distribution'].unique())\n",
    "transp_cost = transp_cost.set_index(['production','distribution']).squeeze()\n",
    "\n",
    "max_prod = pd.Series([180,200,140,80,180], index = production, name = \"max_production\")\n",
    "n_demand = pd.Series([89,95,121,101,116,181], index = distribution, name = \"demand\")\n",
    "\n",
    "# min overall production requirement, and min shipment amount from p to d\n",
    "frac = 0.75\n",
    "C = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code for the 'gurobipy' model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gp.Model('widgets')\n",
    "m.setParam('OutputFlag',0)\n",
    "m.setParam('PoolSearchMode',2)   \n",
    "m.setParam('PoolSolutions',2)   \n",
    "\n",
    "# decision vars\n",
    "x = m.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m.addConstrs((gp.quicksum(x[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m.addConstrs((x.sum('*', d) >= n_demand[d] for d in distribution), 'meet_demand')\n",
    "\n",
    "# objective\n",
    "total_cost = gp.quicksum(transp_cost[p,d]*x[p,d] for p in production for d in distribution)\n",
    "m.setObjective(total_cost, GRB.MINIMIZE)\n",
    "\n",
    "# optimize and extract and show solution\n",
    "m.optimize()\n",
    "\n",
    "x_values = pd.Series(m.getAttr('X', x), name = \"shipment\", index = transp_cost.index)\n",
    "sol = pd.concat([transp_cost, x_values], axis=1)\n",
    "\n",
    "m_ObjVal = m.ObjVal\n",
    "print(f\"This model has a total cost of {round(m_ObjVal,2)}\")\n",
    "sol[sol.shipment > 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Piecewise Linear Modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success alert-dismissible\">\n",
    "  <a href=\"#\" class=\"close\" data-dismiss=\"alert\" aria-label=\"close\">&times;</a>\n",
    "    <strong>Here are a couple a good documentation links to help you use PWL in the code below. </strong>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Links to PWL in `gurobipy`:\n",
    "- [Model.addGenConstrPWL](https://docs.gurobi.com/projects/optimizer/en/current/reference/python/model.html#Model.addGenConstrPWL) command\n",
    "- [Visualization of PWL](https://docs.gurobi.com/projects/optimizer/en/current/concepts/modeling/objectives.html#piecewise-linear-objectives)\n",
    "\n",
    "### Scenario 1\n",
    "We have a new cost structure for transporting the widgets we make. \n",
    "- There will be a flat 0.50 increase for all transportation costs. \n",
    "- If we decide to bulk ship at least half of a production facility's max to any one distribution location, every widget over that half capacity number costs 60% of the new transportation cost. \n",
    "\n",
    " So, if the original cost is 2, the new cost of shipping from $p$ to $d$ for the first half of capacity is 2.5 (0.50 flat increase), and anything produced after that ships at 1.5. Update the model to reflect this (0.60 * 2.5). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_pwl = gp.Model('widgets')\n",
    "m_pwl.setParam('OutputFlag', 0)\n",
    "\n",
    "x_pwl = m_pwl.addVars(production, distribution, vtype=GRB.SEMICONT, lb = 30, name = 'prod_ship')\n",
    "pwl_cost = m_pwl.addVars(production, distribution, name = 'pwl_cost')\n",
    "\n",
    "# constraints\n",
    "can_produce = m_pwl.addConstrs((gp.quicksum(x_pwl[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m_pwl.addConstrs((gp.quicksum(x_pwl[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m_pwl.addConstrs((x_pwl.sum('*', d) >= n_demand[d] for d in distribution), name = 'demand')\n",
    "\n",
    "# create empty series for analysis \n",
    "bulk_costs = pd.Series(index=transp_cost.keys(), name='bulk_costs')\n",
    "base_costs = pd.Series(index=transp_cost.keys(), name='base_costs')\n",
    "\n",
    "# Create piecewise cost functions\n",
    "for p in production:\n",
    "    for d in distribution:\n",
    "        max_produce = ???\n",
    "        frac_produce = ???\n",
    "        \n",
    "        # Adjusted costs\n",
    "        base_cost = ???\n",
    "        bulk_cost = round(???, 2) # 25% increase for the first half\n",
    "\n",
    "        base_costs[p,d] = base_cost\n",
    "        bulk_costs[p,d] = bulk_cost\n",
    "\n",
    "        # Define the piecewise-linear function points\n",
    "        pwl_points_x = [\n",
    "            0, ???, ???\n",
    "        ]\n",
    "\n",
    "        pwl_points_y = [\n",
    "            0, ???, ???\n",
    "        ]\n",
    "\n",
    "        m_pwl.addGenConstrPWL(???, ???, ???, ???, name=f\"pwl_{p}_{d}\")\n",
    "\n",
    "# Objective: minimize the total cost\n",
    "m_pwl.setObjective(gp.quicksum(pwl_cost[p, d] for p, d in transp_cost.keys()), GRB.MINIMIZE)\n",
    "\n",
    "# Optimize\n",
    "m_pwl.optimize()\n",
    "\n",
    "x_values_pwl = pd.Series(m_pwl.getAttr('X', x_pwl), name = \"shipment_pwl\", index = transp_cost.index)\n",
    "sol_pwl = pd.concat([transp_cost, base_costs, bulk_costs, x_values_pwl], axis=1)\n",
    "\n",
    "print(f\"This model has a total cost of {round(m_pwl.ObjVal,2)}\")\n",
    "sol_pwl[sol_pwl.shipment_pwl > 0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_comp = pd.concat([sol, sol_pwl.drop('cost', axis=1)], axis=1)\n",
    "sol_comp[(sol_comp.shipment > 0) | (sol_comp.shipment_pwl > 0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scenario 2\n",
    "The original model does not reflect any production costs to make a widget in different locations. Suppose we have a given production cost per widget for making up to 25% of a production facility's capacity. After that level of production is reached, then the cost *increases* by 25%. For example, if a facility can produce 100 widgets and the initial cost is $1, then the first 25 widgets made will cost $1 each and up to the next 75 widgets qill cost $1.25 each to make. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_pwl2 = gp.Model('widgets')\n",
    "#m_pwl2.setParam('OutputFlag',0)\n",
    "\n",
    "# decision vars\n",
    "x_pwl2 = m_pwl2.addVars(production, distribution, vtype=GRB.SEMICONT, lb = C, name = 'prod_ship')\n",
    "\n",
    "# constraints\n",
    "can_produce = m_pwl2.addConstrs((gp.quicksum(x_pwl2[p,d] for d in distribution) <= max_prod[p] for p in production), name = 'can_produce')\n",
    "must_produce = m_pwl2.addConstrs((gp.quicksum(x_pwl2[p,d] for d in distribution) >= frac*max_prod[p] for p in production), name = 'must_produce')\n",
    "meet_demand = m_pwl2.addConstrs(x_pwl2.sum('*', d) >= n_demand[d] for d in distribution)\n",
    "\n",
    "# Total production for each facility\n",
    "base_prod_cost = pd.Series([1,1.2,.8,.8,.9], index = production, name = \"max_production\")\n",
    "\n",
    "# Add new variables and/or constraints\n",
    "q = m_pwl2.addVars(production, name='q')\n",
    "pwl_prod = m_pwl2.addVars(production, name='pwl_prod')\n",
    "m_pwl2.addConstrs((q[p] == gp.quicksum(x_pwl2[p, d] for d in distribution) for p in production), name='total_production')\n",
    "\n",
    "for p in production:\n",
    "    m_pwl2.addGenConstrPWL(\n",
    "        ???,\n",
    "        ???,\n",
    "        ???,\n",
    "        ???,\n",
    "        name=f'pwl_cost_{p}'\n",
    "    )\n",
    "\n",
    "# objective\n",
    "total_transp_cost = gp.quicksum(??? for p in production for d in distribution)\n",
    "m_pwl2.setObjective(total_transp_cost + ???, GRB.MINIMIZE)\n",
    "\n",
    "# optimize and extract and show solution\n",
    "m_pwl2.optimize()\n",
    "\n",
    "x_values_prod = pd.Series(m_pwl2.getAttr('X', x_pwl2), name = \"shipment\", index = transp_cost.index)\n",
    "sol_prod = pd.concat([transp_cost, x_values_prod], axis=1)\n",
    "print(f\"This model has a total cost of {round(m_pwl2.ObjVal,2)}\")\n",
    "sol_prod[sol_prod.shipment > 0].drop('cost',axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n",
    "\n",
    "The first practical method discussed was basic Monte Carlo simulation. In the code above, we used [solution pools](https://docs.gurobi.com/projects/optimizer/en/current/features/solutionpool.html) to get the top two solutions. We already extracted the optimal solution. The code below queries the second best solution we asked Gurobi to find and puts the solutions together for comparison. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.Params.SolutionNumber = 1\n",
    "m_ObjVal2 = m.PoolObjVal\n",
    "print(f\"This model has a total cost of {round(m_ObjVal2,2)}. The original cost was {round(m_ObjVal,2)}\")\n",
    "x_values2 = pd.Series(m.getAttr('Xn', x), name = \"shipment2\", index = transp_cost.index)\n",
    "sol2 = pd.concat([transp_cost, x_values2], axis=1)\n",
    "sol_sim_comp = pd.concat([sol, sol2.drop('cost', axis = 1)], axis=1)\n",
    "sol_sim_comp[(sol_sim_comp.shipment > 0) | (sol_sim_comp.shipment2 > 0)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in the blanks in the code below to simulate an uncertain transportation cost from Cleveland to Indianapolis. We will assume the cost is normally distributed with mean 2.3 and standard deviation 0.2. Generate 1000 samples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2579)\n",
    "\n",
    "# Number of samples to generate\n",
    "n_samples = 1000\n",
    "\n",
    "# Current cost between Cleveland and Indianapolis\n",
    "mean_cost = 2.3\n",
    "std_sev = 0.2\n",
    "# Generate random costs using normal distribution\n",
    "normal_costs = np.random.normal(loc=???, scale=???, size=???)\n",
    "\n",
    "# Extract the shipment amount from Cleveland to Indianapolis for each solution\n",
    "shipment_amount = sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'shipment']\n",
    "shipment_amount2 = sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'shipment2']\n",
    "\n",
    "# Find the total costs for all other shipments for each of the two solution since we are assuming they are fixed\n",
    "obj1 = sum(sol_sim_comp.cost * sol_sim_comp.shipment) - shipment_amount*sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'cost']\n",
    "obj2 = sum(sol_sim_comp.cost * sol_sim_comp.shipment2) - shipment_amount2*sol_sim_comp.loc[('Cleveland', 'Indianapolis'), 'cost']\n",
    "\n",
    "# Create a DataFrame to hold total costs including the new simulated costs \n",
    "total_costs_df = pd.DataFrame({\n",
    "    'Total_Cost_Sol1': ???\n",
    "    'Total_Cost_Sol2': ???\n",
    "})\n",
    "\n",
    "# Determine the common range for both datasets\n",
    "min_cost = min(total_costs_df['Total_Cost_Sol1'].min(), total_costs_df['Total_Cost_Sol2'].min())\n",
    "max_cost = max(total_costs_df['Total_Cost_Sol1'].max(), total_costs_df['Total_Cost_Sol2'].max())\n",
    "\n",
    "### This is code to create a visualization of each solution's variability\n",
    "# Define bin width and create bins\n",
    "binwidth = (max_cost - min_cost) / 30  # Adjust the divisor for desired number of bins\n",
    "bins = np.arange(min_cost, max_cost + binwidth, binwidth)\n",
    "\n",
    "# Plot the overlaid histograms with the same binwidth and range\n",
    "plt.figure(figsize=(10, 6))\n",
    "\n",
    "# Histogram for total cost using normal distribution costs\n",
    "sns.histplot(total_costs_df['Total_Cost_Sol1'], bins=bins, color='blue', label='Optimal Solution',\n",
    "             alpha=0.5, kde=False, stat='frequency')\n",
    "\n",
    "# Histogram for total cost using triangular distribution costs\n",
    "sns.histplot(total_costs_df['Total_Cost_Sol2'], bins=bins, color='orange', label='Second Best Solution',\n",
    "             alpha=0.5, kde=False, stat='frequency')\n",
    "\n",
    "# Add vertical lines at the means\n",
    "plt.axvline(total_costs_df['Total_Cost_Sol1'].mean(), color='blue', linestyle='dashed', linewidth=2)\n",
    "plt.axvline(total_costs_df['Total_Cost_Sol2'].mean(), color='orange', linestyle='dashed', linewidth=2)\n",
    "\n",
    "plt.title('Overlaid Histograms')\n",
    "plt.xlabel('Cle-Ind Cost')\n",
    "plt.ylabel('Frequency')\n",
    "plt.figtext(0.5, 0.01, 'Dashed lines show respective means', ha='center', fontsize=10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at this histogram, how would you argue going forward with either solution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Individual chance constraint\n",
    "\n",
    "Meeting widget demand in Nashville is now an extremely high priority. Historically, our data says that the demand here for the next sales period follows a [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution) that has a mean of 100 (in the original model, demand was 101). Fill in the code blanks below to make sure we meet the demand with probability 0.95."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import poisson\n",
    "\n",
    "q95 = poisson.ppf(???, ???)\n",
    "q95"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the demand constraint's attribute 'RHS' to quickly update the demand and use the code below to resolve and compare. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "meet_demand['Nashville'].RHS = ???\n",
    "m.update()\n",
    "m.optimize()\n",
    "\n",
    "x_values_Nash = pd.Series(m.getAttr('X', x), name = 'shipment_Nash', index = transp_cost.index)\n",
    "\n",
    "m_Nash_ObjVal = m.ObjVal\n",
    "print(f\"This model has a total cost of {round(m_Nash_ObjVal,2)}. The original model has a cost of {round(m_ObjVal,2)}\")\n",
    "sol_Nash = pd.concat([transp_cost, x_values_Nash], axis=1)\n",
    "sol_Nash_comp = pd.concat([sol, sol_Nash.drop('cost', axis = 1)], axis=1)\n",
    "sol_Nash_comp[(sol_Nash_comp.shipment > 0) | (sol_Nash_comp.shipment_Nash > 0)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scenario Optimization\n",
    "\n",
    "Here we build on the example from the earlier presentation where we were considering three scenarios in portfolio optimization. We'll combine the concept of scenario optimization with \"maxi-min\" from Opti 201. We are considering investing in seven possible sectors under three different scenarios. Below is the input data which defines the sectors, scenarios, and the respective estimated return. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define sectors\n",
    "sectors = [\n",
    "    'Tech',\n",
    "    'TradRetailUtil',\n",
    "    'Renewable',\n",
    "    'FossilFuel',\n",
    "    'Auto',\n",
    "    'GenHealth',\n",
    "    'TradHealth'\n",
    "]\n",
    "\n",
    "# Define scenarios\n",
    "scenarios = [1, 2, 3]\n",
    "\n",
    "# Expected returns (%) for each sector under each scenario\n",
    "estimated_returns = {\n",
    "    ('Tech', 1): 12,\n",
    "    ('Tech', 2): 1,\n",
    "    ('Tech', 3): 3,\n",
    "    ('TradRetailUtil', 1): 3,\n",
    "    ('TradRetailUtil', 2): 4,\n",
    "    ('TradRetailUtil', 3): 3,\n",
    "    ('Renewable', 1): 5,\n",
    "    ('Renewable', 2): 15,\n",
    "    ('Renewable', 3): 4,\n",
    "    ('FossilFuel', 1): 3,\n",
    "    ('FossilFuel', 2): -5,\n",
    "    ('FossilFuel', 3): 4,\n",
    "    ('Auto', 1): 4,\n",
    "    ('Auto', 2): 2,\n",
    "    ('Auto', 3): 4,\n",
    "    ('GenHealth', 1): 1,\n",
    "    ('GenHealth', 2): 4,\n",
    "    ('GenHealth', 3): 10,\n",
    "    ('TradHealth', 1): 3,\n",
    "    ('TradHealth', 2): 2,\n",
    "    ('TradHealth', 3): -3\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an optimization model where we decide the percentage of our budget to invest to **maximize the minimum** total estimated returns considering all scenarios. **HINT**: Define a single new auxiliary variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new model\n",
    "???\n",
    "\n",
    "# Decision variables: Investment fractions in each sector\n",
    "???\n",
    "\n",
    "# Auxiliary variable: Minimum portfolio return\n",
    "???\n",
    "\n",
    "# Add investment allocation constraint\n",
    "???\n",
    "\n",
    "# Add portfolio return constraints for each scenario\n",
    "???\n",
    "\n",
    "# Set objective: Maximize the minimum portfolio return\n",
    "???\n",
    "\n",
    "# Optimize the model\n",
    "???\n",
    "\n",
    "print(f\"Optimal objective value: {round(???,2)}\")\n",
    "\n",
    "# Get the results\n",
    "???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You now decide that your investments need to be a little more diversified. Add any necessary variables and constraints that will:\n",
    "- Make sure you invest in at least four sectors.\n",
    "- If an investment is made in a sector, then it will need to be at least 10% of the budget."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add new variables and/or constraints\n",
    "\n",
    "\n",
    "# Optimize the model\n",
    "\n",
    "\n",
    "# Get the results\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gurobi_ml",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
